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Abstract. - Ubiquitous in nature, convection cells are a clear signature of systems out-of- 
equilibrium. Typically, they are driven by external forces, like gravity (in combination with 
temperature gradients) or shear. In this article, we show the existence of such cells in possibly 
the simplest system, one that involves only a temperature gradient. In particular, we consider 
an Ising lattice gas on a square lattice, in contact with two thermal reservoirs, one at infinite 
temperature and another at T. When this system settles into a non- equilibrium stationary state, 
many interesting phenomena exist. One of these is the emergence of convection cells, driven by 
spontaneous symmetry breaking when T is set below the critical temperature. 



Introduction. — In nature, convection cells can be 
found everywhere, existing at length scales from the mi- 
croscopic to the global. Rayleigh-Benard cells |T] are likely 
the most well known examples. Another are the vor- 
tices (in, e.g., cloud formations) created by the Kelvin- 
Helmholtz instability [2j[3]. In all cases, external forces 
drive convection - temperature gradients and gravity in 
the former case and shear in the latter. Clearly, none of 
these cells can persist under conditions of thermal equi- 
librium. Instead, they are distinct signatures of systems 
far from equilibrium. Motivated by our interest in un- 
derstanding fundamental issues in non-equilibrium statis- 
tical mechanics, we ask here: What are the minimal con- 
ditions for convection cells to persist? The answer may 
provide valuable insight into the essential ingredients as- 
sociated with physics far from equilibrium. To identify 
these conditions, we construct minimal systems, seeking 
the equivalent of the Lenz-Ising model [4] in the study 
of phase transitions. In this spirit, we focus specifically 
on non-equilibrium steady states (NESS), in which con- 
vection cells persist in a time-independent manner. Here, 
we present possibly the simplest model with such proper- 
ties. Unlike the examples above, our system is taken out 
of equilibrium by a temperature gradient alone. In par- 
ticular, we couple two (spatial) sectors of an Ising lattice 
gas to two different thermal baths, evolving with famil- 
iar dynamics and rates. Thus, the underlying dynamics 
fully respect the Ising symmetry. With no other exter- 
nal forces, our cells are the consequence of spontaneous 
symmetry breaking, i.e., phase segregation. 



The properties of the Lenz-Ising model, in its lat- 
tice gas version [5, 6 and in thermal equilibrium, are 
well established. Especially prominent is a second or- 
der phase transition at the Onsager [7J[5] temperature, 
To = 0.5673 J/ks, in half-filled systems in two dimen- 
sions (2D). To drive this model into a NESS and explore 
its behavior far from equilibrium, many authors have in- 
troduced a variety of changes to its usual dynamics, e.g., 
allowing a uniform bias along one of the lattice direc- 
tions [HJUn] or using two different T's for exchanges along 
the two axes [TTM14) . A wide range of novel phenomena 
emerged, leading to a better understanding of NESS in 
general [TSJUB] . Our approach here follows similar lines, 
in that we couple the system to two thermal baths, but in 
a manner that is far more common in nature. As in ear- 
lier studies of NESS, the results are novel and surprising, 
with the most notable being the emergence of convection 
cells. In the remainder of this letter, we define our model, 
present simulation results and quantitative measures for 
these cells. Next, we provide an exactly solvable (small) 
system which displays such cells, as well as steps toward a 
field theoretic framework for analyzing the phenomenon. 
Deferring systematic details to a future publication, we 
close with a discussion and an outlook for future studies. 

Model specifications. — Our system consists of a 2D 
Ising model with the usual nearest-neighbor (NN) inter- 
action on a L x x L y square lattice (with L x and L y both 
even, for simplicity). In a lattice gas, a site (x, y) may 
be occupied by a particle or left empty, i.e., n(x,y) = 1 
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or 0. The interparticle interactions are attractive, i.e., 
H = —J^2n(x,y)n(x',y') with the sum over NN sites 
and J > 0. The system evolves with ordinary Kawasaki 
dynamics [IT]: A NN pair (also referred to as a "bond") 
is chosen at random and exchanged according to the fa- 
miliar Metropolis rates: min[l, e~ AW / fcBT ], where AH is 
the change in T-i due to the exchange. Specifying the 
boundary conditions and starting with an initial config- 
uration, this system will settle into thermal equilibrium, 
with well-known Ising lattice gas properties. In particu- 
lar, for T < To, a half- filled system displays co-existence 
of phases, i.e., a high-density domain separated from a 
low-density one by well defined interfaces. With periodic 
boundary conditions (PBC), each domain is a strip, gen- 
erally aligned along an axis so that its interfaces span the 
smaller of L x , L y . 

To impose coupling to two thermal baths, we partition 
the lattice into two sectors (say, the "left" and "right" 
half-planes: x < L x /2 and > L x /2) and use different T"s 
for each0 This form of drive, with a macroscopic tem- 
perature gradient, is obviously ubiquitous, from stove-top 
cooking to large scale weather patterns. Specifically, we 
set the right bath temperature, T", to infinity (for sim- 
plicity), but leave the left bath at a controllable T. In 
other words, if a particle-hole pair lying strictly in the left 
sector is chosen, they are updated precisely as in the ordi- 
nary Ising case. Otherwise, the pair is always exchanged 
^ e -A-H/oc _ ^ p an , straddling the sectors, (L x /2,y)- 
(1 + L x /2,y), is also exchanged (again, for simplicity). As 
for the boundary conditions, the simplest would be peri- 
odic (PBC). However, to facilitate the measurement of 
convection cells, we impose pinned conditions at the x- 
boundaries. Specifically, we introduce an extra column 
of immobile objects at x = 0: particles at the "bottom" 
(ye [I, Ly/2]) and holes at the "top" (ye [1 + L y /2, £„]). 
The y-boundaries are periodic. Below, we will discuss how 
to detect these cells in a system with full PBC. Here, let 
us emphasize that the pin introduced here plays an iden- 
tical role for an Ising model in equilibrium, in which the 
average magentisation of a finite system with full PBC is 
zero for all T. 

In our simulations, a Monte Carlo Step (MCS) consists 
of L x L y attempts to exchange pairs. Starting with a ran- 
dom, half- filled configuration, we discard 1.6 x 10 6 MCS 
so that the system has relaxed into a steady state. The 
next K MCS are designated as our "run" . Mostly, we have 
K = 3 x 10 7 and carry out 80 runs for better statistics. In 
addition to the usual measurements of the average density 
(n (x, y)), we record every successful exchange over the en- 
tire run. The latter allows us to define net currents across 
each bond, e.g., j x (x, y) is the total number of traverses by 
a particle from (x, y) to (x + 1, y), minus such traverses by 
a hole, divided by K. As a check, for each site we compute 
div j = j x (x, y) -j x (x - 1, y)+i y (x, y) -j y (x, y - 1) and 

1 Other non-equilibrium Ising models coupled to two temperature 
baths, especially using Glauber spin flip dynamics, have also been 
studied. For a review and references, see e.g., Section VII-B in |15| . 



see that it is precisely [n (x, y; K) — n (x, y; 0)] / K . Thus, 
it vanishes in the K — > oo limit. 

By contrast, there is no such constraint on the vorticity 
ijj = curl j (a scalar in 2D). On the lattice, u> (x, y) can be 
associated with the plaquette centered at (x + h, y + | ) . 
A discrete version of § j ■ d£ around the square, 

w (x, y) = j x (x, y) + j y (x + l,y)- j x (x, y + 1)- jy («, V) 

(1) 

is positive for counter-clockwise circulation. For an equi- 
librium system, the average current, (j), must vanish, and 
so, (w) = 0. Being stochastic, uj(x,y;t) should be a ran- 
dom walk in time. Thus, with our definition of j, (uj 2 ) is 
expected to vanish as 1/K. For a system in a NESS, how- 
ever, probability currents are non-trivial in general |18lll9j . 
so that both (j) and (uj) may also be non-trivial. Of 
course, in our model, we do not expect the presence of 
uniform, global currents or vorticity (such as in the driven 
lattice gas [HHUI]). Instead, the symmetries dictate the 
presence of convection cells in equal and opposite pairs. 

Simulation results. — Focusing on the existence of 
oj, we defer our systematic study to a future publica- 
tion and present only one set of results here, namely, 
T = 0.88T o with L x = 20,50,100 and L y = 
20, 50, 100, 200, 400. Fig. la shows a typical configuration 
in a 100 x 100 system and Fig. lb shows a plot of (n (x, y)). 
Not surprisingly, we find phase segregation in the left sec- 
tor, with a high-density domain "at the bottom" (due to 
the pinned BC). Also expected is a homogeneous, half- 
filled state on the right. Across the sectors, there are den- 
sity gradients which tend to drive particles/holes from the 
high/low density domain on the left to the right sector. 
This leads to emerging local excesses in the right sector 
which, in turn, are also unsustainable, due to the under- 
lying dynamics, so that we expect nontrivial j, uj, and 
persistent convection cells. Such a scenario is sketched in 
Fig lc. 

To quantify j and uj, we first verify that the net number 
of exchanges across each bond indeed increases linearly 
with K for each run. Averaging over the 80 runs, we es- 
tablish a non- vanishing (j (x, y)). To display a vector field 
is typically cumbersome. However, in 2D, a convenient 
quantity is the scalar stream function ip(x,y), the curl 
of which is j. In Fig. 2a, we show if> for a 50 x 50 sys- 
tem, illustrating that the currents are highly nonlocal, i.e., 
spread out over the entire lattic^! For easier interpreta- 
tion, we note that j x = d y rp and j y = —d x ip. Hence, j 
can be visualized from -0 by taking its gradient and ro- 
tating the resultant vectors by 90°, so that the peak at 
x = y = 25 corresponds to the center of the uj > vortex 
(red online) in Fig lc. Another perspective is provided 

2 In order to appreciate this remarkable observation, let us em- 
phasize that the currents would be identically zero everywhere, had 
the two temperatures been the same. Our dynamics violates de- 
tailed balance at a local level (three columns of exchanges), yet the 
currents - a prominent signature of a non-equilibrium system - are 
nonlocal. 
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Fig. 1: (a) Typical configuration of a 100 x 100 system, with 
T £S 0.88Tb and T' = oo. (b) The average density profile, 
(c) A sketch of a 20 x 20 case, illustrating the presence of the 
convection cells. For clarity, this system is displaced "upwards" 
by Ly/4. Thus, the top row shown here is the y — 15 line. 



by the associated u, shown in Fig 2b. Notably, it is es- 
sentially localized at the sector boundary (x — 25), with 
peaks/dips centered on the convection cells. Within each 
sector, just one or two columns away from the boundary, 
uj drops precipitously to noisy backgrounds. To highlight 
this localization and its secondary characteristics, we illus- 
trate w in a 20 x 400 system (sector boundary at x = 10). 
In Fig. 2c, we plot u (x, y) against y for several values of 
x. Note that (i) the width of the primary vortex in y is 
only about 10, (ii) the width in x is O (1), and (hi) there 
are two secondary, counterfiow cells here. 

The contrast between tp and U) should caution us on 
how best to describe "non-equilibrium effects" - some as- 
pects being system-wide and others appearing to be lo- 
calized. Of course, such constrasts are well known: In 
electrostatics, the potential of a localized charge ranges 
throughout space. Indeed, the analogy here is precise, 
since V 2 ip = — w. To emphasize the absence of vortices 
in an equilibrium system, we consider a model which is 
superficially similar to ours, by setting J — on all bonds 
centered on x > Lx 2 +1 , but updating with a single T every- 



Fig. 2: Persistent currents and vortices present in a 50 x 50 
system, with T = 0.88Tb and T' = oo, characterized by (a) the 
stream function ip and (b) the vorticity u>. (c) Vorticity in a 
20 x 400 system: u) (x,y) vs. y for several values of x, showing 
that it is essentially localized at the sector boundary (x — 10), 
with counter-flow cells about 10 lattice spacings on either side 
of the primary vortex. Arbitrary units are used for both ifr and 
w. 

where. The dynamical rules within the right sector here 
are identical to our model above. The only difference be- 
tween the two lies with the updates on the three sets of 
bonds centered on i = y, Lj: 2 ±1 . Simulations (to be re- 
ported elsewhere) indeed show all w's distributed around 
zero, with a width vanishing as K^ 1 ^ 2 . 

Theoretical considerations. — To formulate a mi- 
croscopic theory for our model is facile: A master equation 
can be written in a line or two. However, to find the sta- 
tionary distribution is far from trivial. Nevertheless, we 
can offer some insight by considering the exact solution 
for the simplest possible case, namely, a 2 x 2 lattice with 
just 2 particles. Now, for just two rows, we impose open 
BC in y (since PBC would be meaningless) and so, just 
one bond is relevant, for coupling and exchange across the 
rows. With pin BC in x, particle-hole exchanges can take 
place on only 4 bonds. Thus, their associated (j)'s in the 
steady state must share the same magnitude and only a 
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single u) exists. With just 6 configurations (i = 1, ...6), it 
is straightforward to write the master equation and find 
the stationary probabilities exactly. From here, we can 
compute the frequency of any exchange to take place and 
arrive at (u) = 2 (l - e - iJ l k ^ T ) / (7 + be~ iJ / kBT ). To 
appreciate the significance of having two baths, we repeat 
this computation with T' < oo and find that the first term 
in the numerator is e _4J / fesT . Hence, (w) vanishes when 
T' = T, as expected. 

A concern may be raised that, together with the open 
BC, the pinned column affects particles much like grav- 
ity and so, our convection cell is no different from those 
observed daily in the atmosphere. Indeed, for the 2x2 
system, the pin, along with the extremely small size, has 
serious consequences, such as the vortex persisting for all 
T < oo. To untangle the effects of pin vs. size, we turn 
to slightly larger systems, of size 2 x I with t < 6, and 
open BC at both (left and right) edges. Now, we can 
re-impose PBC in y so that translational invariance is 
restored. Of course, this delocalizes the vortex and the 
average vorticity, (cj), is necessarily zero. Nevertheless, 
despite the absence of a gravity-like pin, convection cells 
are still present and can be detected through a density- 
vorticity correlation, such as (nw). Using exact numerical 
techniques [20], non-trivial values for in (1, y) u> (1, y)) and 
(n (1, y) lj (1, y + 1)} emerge, persisting again for all finite 
T. We believe that this behavior has a counterpart in an 
equilibrium Ising model, e.g., in the two-spin correlation 

(s (x, y) s (x + Tp, y + -rf^ )■ This quantity is non-zero for 
all T < oo for a finite system. The difference between the 
two sectors T < To and T >To emerges only for L — > oo: 
the correlation vanishes in the latter and remains non-zero 
for T < To- Of course, a systematic finite-size analysis of 
(nuj) will address these issues decisively. Simulation stud- 
ies for this purpose are in progress. 

Small systems, though exactly solvable, are clearly of 
limited value for understanding large scale collective be- 
havior. For the latter, a standard theoretical approach 
is to exploit Langevin equations for a mesoscopic con- 
tinuum particle density p (x, t) . For an Ising lattice gas 
evolving towards thermal equilibrium, a well established 
route is model B [3T]: d t p oc V 2 (SF/Sp) + -q, where T 
is the Landau-Ginzburg free energy functional and i] is a 
conserved noise. In this approach, the vorticity never ap- 
pears, since the right hand side is just the divergence of 
the current. Although curl j can never affect the evolution 
of p, it can be significant, especially for non-equilibrium 
systems. However, if we take the form for j from model 
B above, then (the deterministic part of) curl j vanishes 
identically and cannot be captured. Instead, it is neces- 
sary to include the mobility factor, a, in the functional 
form of j : 

JM=^[P]{-V^j . (2) 

Such a factor has been exploited in other studies, espe- 
cially in systems with external drives like gravitational 



or electric fields [33] [23]. Typically, a is assumed to be 
p(l — p) for the lattice gas [THOU]. To have curl j ' ^ 0, 
a further ingredient is necessary, namely, explicit spatial 
dependence in either a or J 7 , or both. Otherwise, curl 
j will be proportional to Vp x Vp = 0. In our case, we 
would model the two sectors coupled to baths with T ^ T" 
by having an explicit x-dependence in J- [p. x]. (It is also 
simple to consider models with a [p, x] , which would in- 
duce shear when driven.) One advantage of this approach 
is that we expect VT (x) to appear, so that curl j will be 
localized to the sector boundary (for the model presented 
here). Work is in progress to turn this intuitive picture 
into a developed mesoscopic theory. 

Summary and Outlook. — In this letter, we report 
the presence of convection cells in a minimal model: an 
Ising lattice gas coupled to two thermal baths. Unlike 
cells common in nature, they are not driven by external 
forces such as gravity or shear. Instead, the microscopic 
dynamics obeys the full Ising (up-down, or particle-hole) 
symmetry, so that the convection cells here are induced 
by the spontaneous breaking of this symmetry (and the 
ensuing density gradients). Through simulations and an 
exact analysis of small systems, we investigate a natural 
quantity to describe convection cells - the vorticity, Ui, or 
curl of the current density. To our knowledge, uj has never 
been studied previously in systems where particles are not 
associated with momenta. We believe our findings offer a 
new perspective on the vorticity. 

A range of interesting questions natural arise from this 
study, some of which can be explored readily. Are the 
vortices affected by the details of the exchanges across 
the sector boundary? How does the spread of ui along 
the boundary (O(10) in Fig 2c) depend on L x l Are the 
properties of the interface (between high- and low-density 
domains on the left) similar to those in equilibrium, or 
are they drastically different? Are there any novel critical 
properties associated with cj, which vanishes along with 
the spontaneous symmetry breaking at high T? What is 
the best characterization of the energy flux from the right 
to the left sector? Then there are, of course, deeper issues 
which deserve careful thought and serious effort. What 
is the origin and behavior of the secondary counterflows? 
And is there a series of smaller and smaller vortices? How 
do the phenomena found here behave in the thermody- 
namic limit? Can the density profiles (e.g., in Fig lc) be 
quantitatively understood? 

Apart from vorticity and convection cells, our two- 
temperature model displays many other intriguing and 
counterintuitive properties. In particular, we found an ex- 
tremely surprising phenomenon in a fully periodic, T > To 
system, in which only a strip of a few columns is coupled to 
T' = oo, namely, the particle density within the strip being 
a bi-modal distribution (symmetric around 1/2)! So far, 
this remarkable behavior is yet to be understood. From 
such surprises and the presence of convection cells, we can 
expect this deceptively simple two-temperature lattice gas 
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to present us with further puzzling challenges. We be- 
lieve that solving these puzzles will provide critical clues 
to formulating an overarching theoretical framework for 
non-equilibrium statistical mechanics. 
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